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A modified method is presented to generate artificial magnetic turbulence that is used for test-particle sim- 
ulations. Such turbulent fields are obtained from the superposition of a set of wave modes with random 
polarizations and random directions of propagation. First, it is shown that the new method simultaneously 
fulfils requirements of isotropy, equal mean amplitude and variance for all field components, and vanishing 
divergence. Second, the number of wave modes required for a stochastic particle behavior is investigated by 
using a Lyapunov approach. For the special case of slab turbulence, it is shown that already for 16 wave 
modes the particle behavior agrees with that shown for considerably larger numbers of wave modes. 
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I. INTRODUCTION 

Although it appears to be a simple problem, the ques- 
tion of how energetic charged particles are scattered 
in turbulent electromagnetic fields is still far from be- 
ing solved. Over the years, a variety of analytical 
theories^"— and (semi-) numerical methods^"— has been 
used. Considerable progress has been made toward un- 
derstanding specific scenarios such as cosmic-ray diffu- 
sion in the Solar systemj^i^ While numerical simula- 
tions can relatively easy be extended for other parameter 
regimes, such is considerably more difficult for analyti- 
cal calculations which have to take into account dif- 
ferent energy regimes, different turbulence models, and 
different time dependences of the turbulence. Further- 
more, wave-particle interactions also play a role in the 
generation of strong magnetic turbulence via plasma in- 
stabilities and subsequent particle radiation, ^^■^■^ both in 
the astrophysical context as well as in laboratory plas- 
mas and fusion devices. However, here the focus will be 
set on test particles, which have a sufficiently low den- 
sity so that the fields remain unaffected. Self-consistent 
calculations have been done as well, which however focus 
on qualitative results based on non-linear perturbation 
approaches'^ or on the combination of mostly irreconcil- 
able ansatzesii^ 

As a basic test case, much effort has been focused on 
magnetostatic fields,— which conserve the kinetic en- 
ergy and where only the particles' direction of motion 
is changed through scattering processes. As a numerical 
counterpart to analytical calculations, special attention 
has been diverted to a Monte-Carlo method, which traces 
the trajectories of test particles in a prescribed artificial 
turbulence. In the case of the Padian codcr^ normalized 



variables are introduced as r = f2t, where f2 = qBQ/{mc) 
is the gyrofrequency, and R = jv/ (flio), which is named 
a normalized rigidity with £o the turbulence bend-over 
scale. The equations of motion then read 
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where and ess are unit vectors denoting the direc- 
tions of the mean and turbulent magnetic fields, respec- 
tively. While the background field is usually homoge- 
neous (sometimes, special geometries are assumed such as 
the Parker spirali^i^S or adiabatically focused fields^^:^^), 
the spatial (and temporal, for example due to plasma 
waves?- ) variation of the turbulent fields are consider- 
ably more complicated. 

Two basic onsets have been used over the years: First, 
the grid approach discretizes both the spatial and the 
wavenumber coordinates to enable a fast Fourier trans- 
form. Before integrating particle trajectories, the turbu- 
lent fields are specified on the grid using a superposition 
of (standing) wave modes. The advantage is that the tur- 
bulence structure has to be evaluated only once, while, 
on the negative side, up to 2^^ Fourier modes-^"^^ have 
to be used in order to obtain a sufficient resolution of 
spatial structures. 

Second, following iGiacalone and Jokipiil i^ (which will 
hereafter be denoted as G J99) the continuous approach is 
based on a just-in-time evaluation of the turbulent fields 
at each particle position, which has the advantage that 
no grid structure is needed. For the test case of isotropic 
turbulence, the three basic requirements for the turbu- 
lence generator are: (i) all three turbulent magnetic field 
components have the same magnitude on average; (ii) 
the probability for the random orientation of the sam- 
ple Fourier modes is isotropic; and (iii) the divergence 
of the turbulent magnetic field must be zero. As shown 
recently,^^ this was not possible using previously used 
approaches because at least one condition was always vi- 
olated. Therefore, we feel that it is appropriate to present 
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Table I. An incomplete survey of wave modes and 
simulations used type in the literature 

Reference No. Simulation type Number of wave modes 



36 


Continuous 


18 


37 


Continuous 


24 


38 


Continuous 


50 


39 


Continuous 


64 


40 


Continuous 


100 


41 


Continuous 


128-1024 


42 


Continuous 


384 


43 


Continuous 


768 


44 


Continuous 


5000 


45 


Grid 


7168 


46 


Grid 


8200 


47 J8 


Grid 


~ 10^ 


49 


Grid 


1,657,864 


24-26 


Grid 


221_222 



a modification of the original method that fulfills all three 
conditions. 

A third, relatively novel ansatz is based on the appli- 
cation of stochastic differential equation solvers^®'^^ for 
the cosmic-ray transport equation. The method has been 
applied to pick-up ions, particle acceleration, and cosmic- 
ray diffusion.^.1^"^ 

For the two first approaches, there appear to be no sys- 
tematic investigations of the number of wave modes that 
is necessary in order to have a sufficiently turbulent elec- 
tromagnetic field structure. While, from a basic point of 
view, isotropic and homogeneous turbulence is obtained 
only in the limit of infinitely many wave modes there 
are indications that anisotropic structures might be oblit- 
erated already for a finite number of wave modes In 
the literature, various numbers of wave modes have been 
used with seemingly no systematic motivation besides a 
fondness for even numbers (see Table HI. In the Padian 
code, usually 128-2048 modes are used. 

The purpose of this article is thefore two-fold: In 
Sec. ini a new method for the generation of turbulent 
electromagnetic fields is presented and the properties of 
the resulting diffusion coefficients are discussed for the 
special case of isotropic turbulence both by direct evalu- 
ation of the magnetic field components and by analyzing 
the results of a test-particle simulation. In Sec. IIIIl a 
systematic investigation in terms of a Lyapunov ansatz 
is undertaken in order to determine the number of wave 
modes that is required to obtain a stochastic particle mo- 
tion in slab turbulence. In Sec. IIV[ the results are sum- 
marized and future applications are discussed. 



II. GENERATING TURBULENCE 

In this section, a general method for the generation of 
artificial magnetic turbulence will be presented, which is 
based on the work of GJ99. 



A. General Method 

The general prescription for the numerical generation 
of magnetic turbulence is that of a three-dimensional, 
stochastic Fourier transform 

esB = £^ ^nM^n) COs[fc„z' + Cn - ^(fcn)*] , (2) 
n=l 



where e is a correction factor (see Sec. Ill B| ). Furthermore, 
C„ is the random phase angle, and uj{k) is the dispersion 
relation if plasma-wave propagation effects are to be in- 
cluded (for magnetostatic turbulence, uj can be set to 
zero). According to GJ99, the polarization vector ^„ is 
given by 



/.GJ99 

4 — I « sm a cos 



—i sin a sin (j> + cos a cos 9 cos < 

i -I- cos a cos 6* sin (3) 
sin 6 cos a 



with additionally exp[j . . . ] instead of cos[. . . ] in Eq. ([2|). 
For simplicity we omitted the index n here but bear in 
mind that the angles a, 9, and (f) are different for each 
wavemode. The sum in Eq. ([2]) extends over Nm logarith- 
mically spaced wavenumbers so that Akn/kn is kept con- 
stant. Throughout, normalized wavenumbers are used as 
k = £ok with k the true wavenumber and £q 0.03 a. u. 
the bend-over scale. 

The amplitude function, A(kn), is given through 



'N,y 



A^kn) = G(fc„)Z\fc„ ^(^-) 



(4) 



where G{k) is the turbulence power spectrum. A useful 
form for G{k) that allows for variable energy and incrtial 
range spectral indices is-52. 



Gik) = 



k'J 



2{SBf£or{{s + q)/2) 

r((s-i)/2)r((q + i)/2) + 



(5) 



'^^ F. AIouani-Bibi, private communication (2012). 



where usually s = 5/3 is chosen to reflect the Kolmogorov 
inertial range. ^^'^^ Note that, due to the normalization 
factor in Eq. (|4]), it is sufficient to use only the last 
factor — i.e., the two power-laws — of Eq. ([5|) when using 
the spectrum for the purpose of numerical turbulence 
generation. 

Furthermore, the use of a cosine instead of combining 
a complex exponential function with a complex polar- 
ization vector has the additional advantage that the re- 
quired computation time can be significantly decreased 
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(by almost a factor 2). Such is owed to the fact that, in 
calculating the particle trajectories, typically more than 
90% of the computation time is spent in the evaluation 
of Eq. (21). Therefore, any optimization in the implemen- 
tation immediately pays off. 

The connection with the magnetic turbulence strength 
is obtained from the condition 



(SB)' 



d3fcTrP^,„(fe), 



(6) 



where the trace of the magnetic correlation tensor, 
Pfm(fc), is integrated. The well-known isotropic formii^ 
reads 



G{k) f ktk 



87rfc2 



fc2 



icr{k)eimn ~jr I 1 (J) 



where a is the magnetic helicity. Inserting Eq. ([7]) in 
Eq. (O yields 



3, G{k) 



Ank^ 



G{k) dk 



(8) 



if G depends solely on the magnitude of the wave vec- 
tor. Such underscores the validity of Eqs. dH) and ^ 
without the need for a modification of the spectral index 
depending on the number of dimensions i^i^ 




FIG. 1. (Color online) Illustration of the angles used for the 
construction of the rotation matrix. Values are used as (j> ~ 
60°, e = 40°, and a = 40°. 



vectors from a rotation in the x'-y' plane: 



B. A New Transformation Matrix 

In previous implementations,^* the rotation matrix as 
given through Eq. ^ [see Eq. (5) of GJ99] had been used. 
Here, in contrast, the generation of random directions 
proceeds as follows. First, with k — kK, a random wave 
(unit) vector is chosen as 



^1 — ?f COS(j)^ 



rf' sin ( 

■n 



(9) 



with (f) e [0,27r] and cos6' = r] E [—1,1] the cosine of 
the polar angle 9. The parameter ij is chosen so that the 
density di? = sin9d9d(j) = dijdcj) is constant over the 
entire sphere, which requirement would not be met if 9 
were to be uniformly distributed in the interval [0, tt]. 

The next step is to choose two vectors that are always 
perpendicular to each other and to /t. A straightforward 
choice is to restrict one vector, y' , to the x-y plane (see 
Fig. [J). To complete the (right-handed) orthogonal sys- 
tem, set x' = y' X K so that 

x' ~ (cos 6 cos (f), cos 9 sin (j), — sin 9) 
y' = sin 0, cos (/>, 0) . 

Obviously, the vectors x' and y' are not random as 
y' was restricted to the x-y plane. Therefore, an addi- 
tional random angle a € [0, 27r] is used to obtain two new 



sm sm a -\- Tj cos (p cos a 

— cos (f) sin a + rj sin cj) cos a 

— ^\ — -rf- cos a 

— sin (\) cos a + r] cos sin a ' 
cos 4> cos a + 7] sin (f> sin a 



(10a) 



(10b) 



sma 



This rotation ensures that all transformed unit vectors 
■0, and K (in particular: xp, By ^ ^, and — >■ k) 

are perpendicular to each other and have random direc- 
tions. Note that the results agree with the "y- convention" 
of the Eulerian rotation matrix M(i/j, £,/«); the random 
parameters are (the Eulerian angles) 0, 77, and a. Note 
also that the choice of Eq. (jlOb[) as the designated po- 
larization vector ^ is completely arbitrary. For compu- 
tational/numerical purposes 1/' and ^ are equivalent and, 
thus, interchangeable. In fact, when compared to the 
polarization vector of GJ99, we find that |C!J99 ^ 

(jlOap are almost identical. The only difference is the 
imaginary number in the x- and y-componcnts of ^ 
®- 

This improved transformation matrix is the first main 
result in this paper. The advancement of this matrix 
is that for the test case of isotropic turbulence all three 
basic requirements (see Introduction) are fulfilled simul- 
taneously: 

a. Isotropy Requirement: The turbulence is, in fact, 
isotropic. Thus, all values of the polarization vector ^ 
have the same probability, so that the expectation value 
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FIG. 2. Sampling of 10** magnetic field values at random positions. Shown are the histograms of the SB^, 5By, and SB^ values 
(upper panels) together with the histogram for the (SBx)^, (SBy)'^, and (SBz)'^ values (lower panels). 



(mean) for each vector component vanishes, 

^=^^ daj^dr^l^ cl(/.C. = 0, (11) 

where i S {x, y, z}. 

b. Divergence Requirement: The divergence of the 
turbulent magnetic fields is zero, V • B = 0. Thus, the 
turbulence is divergence free, which corresponds to fc-^ = 
0. This condition is always fulfilled, since fc and ^ were 
constructed that way. 

c. Magnitude Requirement: All three turbulent 
magnetic field components have on average the same 
magnitude, 

''^Bi^l ""//v. '^'-''-r 

which, according to Eq. ([2]), leads to {SBiY = 1/3. 

At this point, we want to emphasize that all require- 
ments are fulfilled simultaneously. This supersedes the 
need for different correction factors, which had to be 
introduced^^ for the classic GJ99 rotation matrix, where 
either ^2 _ ^2 ^2 qj.^ jg corrected, the turbu- 

lent field is not divergence-free or, if that is corrected, 
the surface density for the wave vector directions is not 
uniform. 

The turbulent magnetic field vector as a function of 
space and time is then determined from Eq. ([2|). There, 
in addition to Eq. (fT2|) . the square average of the cosine 
contributes a factor l/\/2, which is remedied via the cor- 
rection factor £ — ^/2. 

C. Numerical Test: Results for Isotropic Turbulence 

As a quick test, we calculate numerically the mean 
squared values of each component of the turbulent mag- 



netic fields and the scattering mean-free paths without 
any guide magnetic field. 

First, numerical sampling of lO'* turbulent magnetic 
field vectors with random spatial positions yields 

{SB^y = 0.342 ± 0.481 
{SBy)^ = 0.314 ±0.443 
(SB,)^ = 0.345 ±0.485. 

The total turbulent magnetic field strength is 

((55)2 ^ 1.00036 ± 7.67 x lO'^ 

in agreement with Eqs. ([2]), dH), and The distribu- 

tion of SBx , SBy , and SB^ together with their squares are 
shown in Fig. [2] It is thereby illustrated that all mag- 
netic field components equally have a mean of zero (up- 
per panels) and approximately equal variances or mean 
square values (lower panels), which corresponds to the 
first and second moments, respectively. Due to the ro- 
tational symmetry inherent in all turbulence generators 
discussed here, it is particularly useful to consider the 
behavior of the z component as opposed to the x and y 
components. 

Second, Fig. [3] shows the results for the "running" 
mean-free paths, 

A,(0 = ^((Z\n)'), ^e{x,y,z}, (13) 

where the angular brackets (• • • ) denote the averaging 
over - in this case ~ an ensemble of 10^ particles. 

As expected, in the long time limit A^, = Ay = Az 
(black lines). Furthermore, three different particle en- 
ergies are assumed via the normalized rigidity, R = 
jv/{iof^)- While for small times, there are some devi- 
ations between the three components, there is clear con- 
vergence towards a common value for large times. 
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FIG. 3. (Color online) The three running mean-free paths 
(solid lines: A^,; dotted Unes: A^; dashed hnes: A^) for the 
rotation matrix described in Sec. |ll] (black lines) as opposed 
to the previously used version^'i* (red lines). The parame- 
ter K = 'yv/[ioQ) denotes the normahzed rigidity, which is 
proportional to the particles' velocity. 

In contrast, the GJ99 rotation matrijs^ii^ exhibits a 
clear difference between the mean-free paths. For all en- 
ergies considered we find that > Xx = \- Somewhat 
interestingly, for almost all particle energies (except high 
energies) the new, corrected, mean free path lies between 
^GJ99 A^y^^. However, for high-energetic particles 
we find that Az was underestimated, while Aa; and Aj^ are 
almost in agreement with the new, common value. 



D. Special Geometries 

The above considerations are valid for the cases of 
isotropic and continuous anisotropic turbulence geome- 
tries, where the latter is for example given by the Maltese 
cross geometry»^i^ However, several additional standard 
cases have to be considered. 

a. Slab Geometry: Slab turbulence is characterized 
by SB{r,t) = SB{z,t), so that k = ke^ and SB _L e^. 
Hence, rj has to be set to unity but otherwise the above 
considerations remain unchanged. Without loss of gen- 
erality we set a = 0, so that the vectors can be rewritten 
as 

^p = (cos (j), sin (j), 0) (14a) 
^ = (-sin</),cos0, 0) (14b) 
/«= (0,0,1). (14c) 

b. 2D Geometry: The opposite case of two- 
dimensional turbulence is characterized via 6B(r,t) — 
SB{x,y,t) so that k _L e^. In addition^i^ one usually 
requires 6B _L e^. Therefore, without loss of generality 



one sets 77 = and a = so that 

1/^^(0,0,-1) (15a) 
^ = (— sin(/),cos0, 0) (15b) 
K (cos 0, sin 0, 0) . (15c) 

If one were to ignore the condition that SB be perpen- 
dicular to the background field, a may also be chosen 
randomly. 

c. Composite Geometry: Slab and 2D geometries 
can be combined to obtain a quasi-three-dimensional 
magnetic field as 

SB 

comp 

(r, t) — SBs\a,h{z, t) + (5i?2D(a;, y, t), (16) 

which is named composite turbulence. Here, care has to 
be taken that the amplitude functions are modified so 
that Asiab(fc) -J> Asiab(fc)?' and A2D{k) A2D(fc)(l - r), 
where r describes the relative weight of the slab contri- 
bution according to {SBY = r{SBs\ahY + (1 — r){5B2ViY- 
Typical values^i^ vary between r — 15% and r = 20%. 

E. Electric Fields 

If plasma waves are to be included via the dispersion 
relation Luik) in Eq. ([2|), turbulent electric fields are in- 
duced via Faraday's law, which (i) modify the diffusion 
coefficients and (ii) lead to stochastic acceleration22i^ via 
momentum diffusion.— i^i^i^ In general, the solution for 
E of the Faraday equation B = {c/uj)k xE is ambiguous. 
However, for basic magnetohydrodynamic (MHD) waves 
such as (fast) magnetosonic waves and Alfven waves, the 
properties of the electric field vector and its amplitude 
are known. 1 In turn, the two MHD wave types will be 
discussed. 

As a side note, there are subtle issues about the use 
of random phase factors in Eq. ^ for the determination 
of particle acceleration. While the assumption of ran- 
domized phases is frequently used in test-particle simula- 
tions as well as for diffusive-shock acceleration, nonlinear 
wave-wave interaction processes tend to generate coher- 
ence among Alfven waveSi®-'-- Therefore, the best way 
of generating the proper electromagnetic field structures 
would be to solve the magnetohydrodynamic equations.— 

1. Fast Magnetosonic Waves 

The relevant properties of the electric field unit vector, 
esE, are: (i) bse -L ess] (ii) esE -L k; and (in) bse is 
independent of Bo- Combining the second condition with 
Faraday's law, the electric field strength is given by 

\E\^^\B\, (17) 

where the dispersion relation w = ±VAk has been used 
with = Bo/y/Anp the Alfven speed. 

Because k x tp = ^, the orientation of the electric field 
is simply that of the vector as given in Eq. (llOap . 
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FIG. 4. (Color online) Trajectories of two particles starting at z = (blue lines) and z = 0.1 io (red lines) together with a 
magnetic field line (black line). From left to right, the three panels show the cases of two, four, and eight wave modes are 
shown, corresponding to the value for the parameter N,n in Eq. ([2|). 



2. Alfven Waves 



A. Simplified Turbulence 



Here, the geometric properties are slightly modified, 
insofar as: (i) i-sE -L ^sb', (ii) ■ k — k&vaO with 
6 = Z{k, Bq) as in Sec. Ill B[ and (iii) esE -L -Bo so that 
the z component is automatically zero. 

Together with Faraday's law, the second condition im- 
phes that, for 9 — 7r/2, the magnetic field decreases and, 
at the same time, the frequency decreases because of the 
dispersion relation uj = zLvAk cosO. Such could be taken 
into account by setting 

^ = (— cos 6* sin (/), cos 6* cos 0, 0) (18a) 
i/j = (cos 0, sin 0, 0) . (18b) 

In that case , the cor rection factor has to be set to e = 
because {SB^^yY = 1/6. However, it is undesirable 
to have a turbulent magnetic field strength that varies 
with the direction and is therefore reduced on average. 
Therefore, it is useful to restrict the wave vector to the 
z direction (i.e., set 6* = 0) so that slab geometry is 
retained. 



III. NUMBER OF WAVEMODES 

A second, related problem is the following: How many 
wave modes (i.e., what value has to be chosen for Nm) 
are needed so that Eq. ([2]) describes a turbulent magnetic 
field in a strict physical sense? In what follows, first a 
criterion is found which is then put to the test by running 
a full simulation. 



To shed some light on this matter, a simplified equation 
of motion, v = {q/c)v x B is solved numerically for a 
magnetic field of the form 



N„ 



B 



{x,y} 



(r,t) = ^^(fc) 



(Q!„)cos(fcz) (19) 



with the amplitude function from Eq. (|4]) and with Q!„ = 
[2{n ~ 1)/N„i + l/4]7r. For a single wave mode, the par- 
ticle motion remains regular. That behavior is already 
modified if Nm = 2. In Fig. |3] the trajectories of two 
neighboring particles for two, four, and eight wave modes 
are illustrated together with a magnetic field line. 

The question whether or not particles are "randomly" 
scattered is investigated by observing the time evolution 
of the distance Az between two initially neighboring par- 
ticles (with a distance of, say, Az — 0.1 Hq). While, Az 
remains constant for one wave mode, an approximately 
linear increase is found for the cases of two and four wave 
modes as exhibited by comparison of the two trajecto- 
ries in Fig. m That behavior is entirely changed for eight 
wave modes, when a stochastic variation is found between 
regimes with positive and negative Az. 

A more rigorous mathematical criterion is found in 
chaotic systems, where initially neighboring instances de- 
viate exponentially according to 



Az{t) sa Z\z(0)exp (At) 



(20a) 
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FIG. 5. The quasi-Lyapunov exponent from Eq. ()20b|) as 
obtained from the time evolution of the spatial deviation be- 
tween two initially neighboring particles. From top to bottom, 
the three panels show the cases of two, four, and eight wave 
modes are shown. Note that, in the upper and middle panels, 
A is negative for early times t < 300. 



where 



A = lini lim — In 

t^oo ziz(0)--i-0 t 



Az{t) 



Az{Q) 



(20b) 



is a quasi-Lyapunov exponent. However, here the total 
particle velocity is constant so that an exponential devia- 
tion would be possible only for small times. Furthermore, 
a truly random behavior leads to diffusion, in which case 
one expects Az cx \/t. Therefore, a quasi-Lyapunov ex- 
ponent should exhibit a behavior A cx ln(i)/i if Nm is 
sufficiently large so that particles diffuse. As shown in 
Fig. [5] for the cases of two, four, and eight wave modes, 
that requirement is almost (but not quite) fulfilled for 
two and four wave modes, while, for eight wave modes, 
an excellent agreement is seen even for small times. The 
fact that A becomes very small at intermittent time in- 




















10"' 








10-=^ 


■lil' = 10-7 








FIG. 6. (Color online) Perpendicular mean-free path, Ax/^o, 
as a function of the normalized simulation time vt/lo- Shown 
are the cases of two (red), four (blue), eight (green), 16 (ma- 
genta dashed), and 512 (black) wave modes. Note that addi- 
tional runs with 32 and 128 wave modes are almost indistin- 
guishable from the case of 512 wave modes. 



tervals underscores that, from time to time, particles can 
move toward each other, as is expected from a random- 
walk motion. 



B. Full Turbulence Simulation 

To verify the above findings, a series of simulations has 
been run using the Padian codei§ for slab turbulence of 
the average strength 6B = Bq. The resulting perpendic- 
ular and parallel mean- free paths are shown in Figs.lHland 
[3 respectively, which have been obtained from averaging 
over lO'^ particles in 250 turbulence realizations. 

The number of wave modes, Nm, as used for the evalu- 
ation of Eq. ([2]) is varied between Nm — 2 and Nm = 512. 
As expected from the previous considerations, diffusive 
behavior is seen for eight wave modes. Both the par- 
allel and the perpendicular running diffusion coefficients 
exhibit the same behavior that is found for a consider- 
ably higher number of wave modes. However, the abso- 
lute value of the mean-free paths still deviate significantly 
from the true values. 

In agreement with the results found in Sec. IIII Al the 
overall conclusion is that eight wave modes are sufficient 
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FIG. 7. (Color online) Same as Fig. (6] only that here the 
parallel mean-free path, A||/^o, is shown. 



for a diffusive particle behavior. The reason is that such 
allows for particle trapping (as exhibited in the right- 
hand panel of Fig.|3]at z = 20 .^o): which however depends 
on the precise initial conditions so that some particles 
become trapped with a subsequent reversal of their par- 
allel motion. Of course, for other than slab turbulence, 
one expects that more wave modes are necessary so that 
the three-dimensional structures are sufficiently resolved. 
Additionally, if different wavenumber ranges contribute 
to the spectrum (such as energy, inertial, and dissipation 
ranges^^), each range has to be sampled with its own set 
of wavenumbers. 



IV. DISCUSSION AND CONCLUSION 

In this article, the generation of numerical electromag- 
netic turbulence has been revisited. Such turbulent fields 
are frequently used in test-particle simulations to deter- 
mine the scattering parameters of energetic particles such 
as cosmic rays. Several methods use the approach of 
superposing a number of randomly propagating Fourier 
modes, which is theoretically well-based.'*^ Whereas pre- 
vious approaches sometimes lacked the required random- 
ness of the Fourier modes as well as the properties of the 
resulting turbulent fields, here a new derivation has been 
presented. Therefore, it has to be pointed out that the 
new derivation presented here solves all three problems 



explained in the introduction:-'^ (i) all three turbulent 
magnetic field components have on average the sample 
magnitude, i.e., {5BiY — 1/3; (ii) if all parameters are 
chosen random, the turbulence is precisely isotropic so 
that a vanishing mean of all turbulent magnetic field com- 
ponents is found, i.e., 5Bi = 0; and (iii) the turbulence 
is divergence-free, which, in Fourier space, corresponds 
to fc • ^ = and is fulfilled by construction. 

In addition, the absolute strength of the turbulent 
magnetic field plays an important role not only for 
the magnitude of the (both running and final, diffu- 
sive) mean free path but also for the entire particle 
behavior as the underlying scaling laws are controlled 
by the Kubo number,— 1^ which is defined 'as K = 
{SB/Bo)/{£\\/i±) if — unlike here — the parallel and per- 
pendicular correlations lengths differ. While previous 
implementations sometimes corrected the transport pa- 
rameters afterwards, the correct advance is to ensure the 
appropriate turbulence strength prior to integrating the 
particle trajectories even if electric fields are included, 
which sometimes reduce the turbulence strength even 
more due to an additional sine factor j^i^ 

Furthermore, the problem of the number of Fourier 
modes required for a sufficiently turbulent magnetic field 
has been tackled from the view point of the particle be- 
havior. It is known that, for slab turbulence, the run- 
ning perpendicular mean free path exhibits a power-law 
as Ax cx while the parallel mean free path is ex- 

pected to be diffusive-SS, Such a behavior is already seen 
for eight Fourier modes [i.e., = 8 in Eq. ([2])]. It 
may seem surprising at first glance that only (figura- 
tively speaking) a handfull of wave modes is sufficient to 
obtain accurate results. But such is in agreement with 
early simulations,—^— which, despite the low number of 
Fourier modes, gave results that were nevertheless ac- 
curate enough to allow for successful comparisons with 
analytical theories. 

In many cases one expects previous simulation re- 
sults to retain their validity due to the robustness of 
the particle transport parameters for certain parameter 
regimes. In general, however, such is not guaranteed for 
the reasons stated above. It is reasonable to assume that 
each additional feature — such as anisotropic correlation 
lengths and/or spectral indices, a dissipation range, or 
time dependence — will require additional modes so that 
the spectrum can be sufficiently resolved. Future work 
should therefore attempt to find criteria as to how many 
wave modes are required depending on the complexity of 
the turbulence power spectrum. 
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